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SECTION  I 


INTRODUCTION 

This  article  studies  some  simple  examples  which  illustrate 
the  working  of  the  Martin-Lomax  iteration  (Ref.  1)  for  the  com- 
putation of  transonic  flow  fields.  They  deal  with  flow  fields 
linearized  for  the  vicinity  of  a fixed  Mach  number  or  (in  the 
discussion  of  the  shock  operator)  linearized  with  respect  to  a 
subsonic  and  to  a supersonic  Mach  number,  both  assumed  to  be  fixed. 
Martin  (Ref.  2)  has  carried  out  similar  discussions  but  only  in 
a rudimentary  form;  this  suffices  to  develop  and  motivate  his 
ideas.  His  discussions  merely  serve  to  suggest  that  a stable 
iteration  scheme  for  flow  computations  can  be  obtained  if  one 
chooses  certain  parameters  in  a suitable  manner,  while  the  actual 
test  of  convergence  is  carried  out  in  a more  realistic  setting. 

For  a justification  of  the  Martin-Lomax  procedure,  the  present 
more  detailed  study  is  therefore  unnecessary.  Nevertheless,  it 
may  have  a certain  interest  because  it  provides  an  insight  into 
the  working  of  the  method.  In  addition,  th^re  are  certain  diffi- 
culties of  an  intuitive  character  which  a novice  to  this  procedure 
is  likely  to  have  and  which  the  author  hopes  to  dispel  by  the 
present  study.  To  be  specific:  in  each  iteration  step  the  Martin 
procedure  solves  a problem  closely  related  to  Poisson's  equation, 
the  solutions  satisfy  the  boundary  conditions  of  an  elliptic 
problem.  It  is  not  immediately  clear  how  the  solution  to  a hyper- 
bolic problem,  which  has  different  boundary  conditions  and 
regions  of  influence  can  be  built  up  in  this  manner.  Of  a 
similar  nature  is  the  question  how  such  a procedure  can  lead  to 
shocks,  for  one  is  inclined  to  ascribe  a smoothing  tendency  to  an 
elliptic  operator,  but  for  a shock  a sudden  change  of  velocity  is 
characteristic.  These  objections  may  appear  less  serious  if  one 
interprets  the  iterations  as  a time  dependent  process,  as  is  done 
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in  Ref.  2.  Then  one  deals  with  a three  dimensional  problem  (two 
space  dimensions  and  the  time)  for  which  the  solution  in  the  long 
time  limit  happens  to  satisfy  the  desired  two  dimensional  prob- 
lems, which  are  governed  by  elliptic  or  hyperbolic  operators. 

Martin's  stability  analysis  of  the  difference  equation  is 
carried  out  by  means  of  a substitute  partial  differential  equation 
which  arises  from  the  three  dimensional  difference  formulation  by 
a Taylor  development  with  respect  to  space  and  time.  Strictly 
speaking,  such  a development  is  valid  only  if  the  wave  length  of  the 
perturbations  is  large  in  comparison  to  the  grid  size,  and  if  the 
changes  with  respect  to  time  are  only  small.  It  is  sufficient  to 
establish  consistency  of  the  difference  scheme  with  the  original 
two  dimensional  partial  differential  equation.  But  the  stability 
analysis  cannot  be  restricted  to  perturbations  of  this  character. 
Short  waves  (down  to  a wave  length  of  twice  the  distance  between 
grid  points)  cannot  be  disregarded.  It  is  true,  in  the  limit  of 
zero  grid  size  such  waves  are  unimportant  in  a smooth  starting 
approximation  for  elliptic  problems.  Still  they  may  be  excited 
by  rounding  errors.  In  supersonic  flows  and  in  flows  with  shocks, 
short  waves  will  be  excited  even  by  the  starting  approximation. 

In  the  author's  opinion  an  independent  analysis  of  the  shock 
point  operator  is  inappropriate.  In  a shock  one  has  a rapid  change 
of  the  gradient  of  the  potential  within  a few  mesh  points.  The 
number  of  mesh  points  where  this  happens  does  not  increase  as  the 
mesh  is  refined.  Moreover,  the  shock  point  operator  is  applied 
only  at  the  shock  points,  that  is  along  a line  of  the  two  dimensional 
flow  field.  It  is  therefore  necessary  that  the  shock  is  considered 
as  imbedded  in  a supersonic-subsonic  flow  field. 

The  present  report  investigates  questions  of  this  kind. 
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SECTION  II 


DESCRIPTION  OF  SAMPLE  PROBLEMS 

Martin  and  Lomax  use  the  velocity  components  as  independent 
variables.  For  the  present  discussion,  it  is  equally  convenient 
to  work  with  the  velocity  potential.  We  consider  the  linear 
differential  equation 


‘ 93T  <Vx>  + ^yy  = 0 (1) 

here  u^  denotes  the  local  deviation  of  the  x-velocity  from  1,  <J>  is  the 
perturbation  potential.  In  our  example  u^  is  considered  as  con- 
stant. In  studying  the  parabolic  operator  it  would  be  necessary 
to  allow  u^  to  change  continuously  with  x.  Equation  (1)  arises, 
if  one  carries  out  a Prandtl-Glauert  coordinate  distortion  for  a 
subsonic  Mach  number  in  a linearized  differential  equatior  for 
subsonic  or  supersonic  flow.  The  local  Mach  number  is  allowed 
to  be  different  from  the  Mach  number  for  which  the  Prandtl-Glauert 
coordinate  distortion  is  carried  out.  The  velocity  fiexd  repre- 
sented by  <j>  gives  a small  perturbation  to  the  local  velocity  u^  . 

The  potential;})  is  assumed  to  be  of  period  L in  the  y direction. 

Let  the  number  of  subdivisions  per  period  in  the  y direction  be 
N,  then  one  has  a grid  size 


h = L/N. 


(2) 


We  have  assumed  that  a Prandtl-Glauert  transformation  has  already 
been  carried  out  in  Equation  (1);  it  is  therefore  reasonable  to 
assume  that  the  grid  size  in  the  x and  y direction  is  the  same. 

In  analyzing  the  procedure  we  consider  particular  solutions 

= (j>U)  -jr)  (3) 
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This  solution  has  n full  waves  within  the  period  L.  In  a finite 
difference  approximation  one  obtains 


/ h 

* - p.  $>Cx)  UhfCn  £-)  4 

- - —7-  /V  *f*(*>{z 

h 


Thus,  one  obtains 


(4) 


Here,  n assumes  values  from  0 to  N/2 ; n = 0 gives  a one  dimensional 
flow.  Let 


^ = (fi  (kk) 


The  difference  form  of  Equation  (4)  for  an  elliptic  problem  (with 
constant  u^  = - | u^  i)  is  then  given  by 

/"//X/  y *r/I  ■ 0 <5> 

Let  the  beginning  and  the  end  in  the  x direction  of  the  region 
under  consideration  be  given  as  k = 0 and  k = K.  The  values  of 
<J>  and  <f>K  are  prescribed  as  boundary  conditions.  For  the  hyperbolic 
region  the  difference  operator  expressing  cf>  xx  is  displaced  by  one 
mesh  in  the  upstream  direction  (here  to  the  left) 

-/%/[ i - ***Yr  («) 


4 

d 


Here  the  values  of  and  x are  prescribed  at  the  initial  point 
(k=0).  To  the  first  order  in  h these  quantities  are  expressed 
by  <j>  0 and<})_1.  Equation  (6)  expresses  6^  in  terms  of 
and  The  initial  conditions  are  of  a nature  which  allows 

the  problem  to  be  solved  by  a marching  procedure. 

A flow  containing  a shock  is  hyperbolic  upstream  (to  the 
left)  of  the  shock  and  elliptic  downstream  of  it.  In  the  Murman 
procedure  (Reference  3)  the  point  at  which  the  shock  operator  is 
to  be  applied  is  found  by  a separate  test.  In  the  present  con- 
text we  deal  with  perturbations  to  a flow  which  already  contains 
a shock.  These  perturbations  are  assumed  to  be  small  enough  so 
that  the  result  of  this  test  will  not  be  affected  by  the  per- 
turbations. Accordingly,  the  shock  point  is  considered  as  fixed. 
Let  s be  the  subscript  of  the  point  at  which  the  shock  point 
operator  is  applied.  To  the  left  and  to  the  right  of  this  point 
one  has  respectively  u£  = |u£|  = const  and  u£  = - | u£ | = const. 

For  1 < k < s - 1 and  for  s + 1 < k < K-l,  one  applies  respectively 
the  hyperbolic  form  (6)  and  the  elliptic  form  (5)  of  the 
difference  operator.  The  shock  point  operator  is  obtained  by 

replacing  3 (u^())x) /3x  by  (h) -1  ( (u£<j>x)  - (u£c(>x)  s_1/2 > • 0ne 

s+1/2 

haS  U£ , s-1/2  = ; US>,s+]/2  = 

and 

($*)  - />''(&*, 

and  because  of  the  shift  by  one  mesh  to  the  left  in  the  hyperbolic 
region 


* h"i -A-j 


One  thus  obtains  for  the  shock  point 
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We  shall  not  carry  out  a study  of  problems  containing  a para- 
bolic point.  In  such  a case  one  would  assume  that  u^  changes 
continusly  from  a negative  to  a positive  value.  Let  k = p be 
the  point  where  the  parabolic  operator  is  to  be  applied.  A simple 
choice  for  u^  is 

It  - Cftti-t  ft-  * ) 

l * 


Then  one  would  have  for  a point  of  the  elliptic  region,  that  is  for 
/ i k i />-/ 

Ce»si  h £ (k  + *£  " p)f\ fat/  -fa)  - fk  - '/*  -/>  Jffa  - 

- if-  Jlx? (H  jj-)  fa  “ 0 

and  in  the  hyperbolic  region,  that  is  for  p + 1 < k < K 

Covii  h {(*,-£  -fa.,)  - fA-*A  -p)(fa_,  -fa-J} 

- * ^ %)fa  “ 0 

At  the  parabolic  point,  one  obtains  by  a procedure  corresponding 
to  that  of  a shock  point 
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e />+  'A 


V*  r i--  - i-  > 


(“fh  ^.,/i  ’%p-i  (f„  - $p-t) 

Hence 

ixtyA)  - 0 

fi 

and  because  of  Equation  (4) 

£/>  - o 

Notice  that  the  parabolic  operator  provides  a boundary  condition 
for  the  elliptic  region.  It  is  therefore  possible  to  compute  the 
elliptic  region  independently.  Of  course,  this  holds  only  for 
the  present  example  where  the  transition  from  the  elliptic  to  the 
hyperbolic  region  occurs  along  a line  x = const. 
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SECTION  III 


THE  MARTIN-LOMAX  ITERATION 


Let  us  denote  by  <)>?;  the  values  <J>.  obtained  by  the  iteration. 

i ^ ^ th 

Assume  that  the  tp,  ' s are  known,  that  is,  that  the  1 iteration  is 

K . i 


completed.  To  obtain  the  equations  for  the  (i+1) 

cima 
i+1 


iteration  the 


. • 2 2 . 
difference  approximation  for  9 4>/9y  is  expressed  in  terms  of  the 


unknown  values  4>, 


while  the  difference  approximations  for  the 


x derivatives  are  computed  for  the 


i. 


s.  Martin  and  Lomax  include 


in  the  equations  further  expressions  which  are  linear  in  the 


evaluated  once  with  the  4> 


i+1 


V 


and  a second  time  with  the  opposite 


If  the  iterations  converge,  that  is,  if 


sign  with  the 

<(>k1+1  = (fi^1  in  the  limit  i -*■  °°,  these  terms  cancel  and  one  is  left 
with  the  difference  approximation  to  the  problem.  These  additional 
expressions  are  of  such  a nature  that  the  resulting  difference 
equations  for  the  <j>k1+^'s  have  constant  coefficients.  Martin's 
preliminary  analysis  suggests  that  for  a suitable  choice  of  certain 
parameters  which  are  left  arbitrary,  the  procedure  will  converge. 
The  equations  ultimately  to  be  solved  are: 


J&r  ejfi'jbfic  poinft 
fir  hyptr l>a//c 
fiy  skoc/i.  f>*i**fo 


(7) 


Similar  equations  arise  in  problems  containing  a parabolic  point 
except  that  then  one  has  to  take  the  fact  into  account  that  u^ 
is  not  constant. 
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SECTION  IV 


DISCUSSION  OF  EQUATION  (7) 

Martin's  discussion  amounts  to  an  application  of  the  von  Neumann 
stability  criterion  to  a simplified  form  of  this  equation.  It  is 
assumed  that  the  region  extends  in  the  x direction  from  - °°  to  + <*> . 
The  elliptic  and  the  hyperbolic  problem  will  be  studied  separately. 
The  difference  equation  for  the  4>^+^'s  has  constant  coefficients. 
Therefore,  it  is  natural  to  study  particular  solutions  of  Eq.  (7) 
of  the  form 


up  (if'*  **  9^ 


(8) 


where  y is  a real  constant,  which  ranges  from  0 to  1/2.  The  smallest 
wave  length  (namely  2h)  occurs  for  y = 1/2.  The  hypothesis  Eq.  (8) 
guarantees  that  these  particular  solutions  remain  bounded  for  k -*■  ±°°. 
The  constant  p^,  so  far  unknown,  is  the  (complex)  factor  by  which 
each  particular  solution  is  multiplied  as  one  proceeds  from  one 
iteration  step  to  the  next  one.  One  obtains  immediately  from  Eq.  (7) 


(f  -ii- 0+1)0- M/jft) 

f/i^  / /-  + 2.  - fcv  cJ/iftt'C  prcl/ems 

I ( f - l Z^vLfh))  jffiy  htptrlo/tc  /trrol/enrs 


From  this  equation,  one  readily  computes  the  p^'s  for  each  choice  of 
0^,  ct2,  lu^l»  n/N,  and  h.  In  order  for  the  procedure  to  converge 
in  a general  case,  it  is  necessary  that  one  finds  values  of  the 
parameters  a-^  and  ct2  in  such  a manner  that  |p^|  < 1 for  all  values  of 
|u^|,  n/N  and  y which  occur  in  the  flow  field. 

To  be  somewhat  more  specific:  one  starts  with  an  assumed  vector 
$k°  which  satisfies  the  boundary  conditions  of  the  elliptic  or 
hyperbolic  problem.  It  can  be  considered  as  a linear  combination  of 
the  desired  solution  (for  which  ^ and  homogeneous  particular 

solution  of  the  form  Eq.  (8).  The  contributions  of  these  homogeneous 


particular  solutions  become  small  as  the  iterations  proceed,  if  |p  | < 1 
for  0 < y < 1/2,  because  they  are  multiplied  in  each  iteration  step  with 
respective  values  of  . Finally,  the  desired  solution  of  the 
inhomogeneous  problem  is  left. 

The  analysis  described  so  far  agrees  with  that  of  Martin  except 
that  he  assumes  n and  y to  be  small.  Such  an  analysis  does  not 
provide  an  insight  into  the  effect  of  boundary  conditions  at  the 
beginning  and  at  the  end  of  a finite  region;  the  conditions  of 
boundedness  at  x -►  ±<»,  which  take  their  place,  are  somewhat  vague. 

For  the  homogeneous  particular  solutions,  one  has  the  boundary 
conditions  <p  = 0 and  <|>v  = 0,  for  hyperbolic  as  well  as  for  elliptic 
problems.  For  the  elliptic  problem,  these  conditions  do  not  offer 
any  conceptual  difficulty,  for  they  agree  with  those  of  the  original 
partial  differential  equation.  For  the  hyperbolic  problem,  one 
postulates  that  <}>q  and  <J>_^  assume  the  values  given  as  initial  condi- 
tions in  all  iterations  (including  the  starting  approximation) ; 
besides  one  sets  $ = 0.  The  boundary  values  for  <p  , and  <j>  are 

used  when  one  computes  the  right  hand  side  for  the  (i  + 1;  iteration 
from  the  values  of  1.  In  the  computation  of  <)>,  1 + ^ only  the  values 

of  $ 1 and  <)>k  appear  as  boundary  conditions. 

If  the  homogeneous  particular  solutions  die  out  in  the  course  of 
the  iterations,  then  the  final  solution  will,  of  course,  have  the 
prescribed  values  of  and  <)>o.  That  is,  it  will  satisfy  the 

boundary  conditions  of  the  hyperbolic  problem.  There  is,  however, 
a question  about  the  meaning  of  the  condition  $ = 0.  Certainly 

there  is  no  room  in  a hyperbolic  problem  to  satisfy  this  condition. 

This  discrepancy  resolves  itself  in  the  following  manner.  It  was 
mentioned  above  that  the  hyperbolic  problem  can  be  solved  by  a 

marching  procedure.  The  expression  Eq.  (6)  allows  one  to  compute 

the  value  of  <f>  for  some  point  k from  points  preceding  k;  it  expresses 
the  requirement  that  the  differential  equation  in  its  finite  difference 
from  be  satisfied  at  point  k.  This  requirement  is  met  for  all 
interior  points  of  the  region,  that  is,  for  l<k<K-l.  At  the 
point  K,  no  such  requirement  is  imposed.  It  follows  that  the  value 
of  4>  = 0 imposed  in  computing  a new  iteration  need  not  coincide 

with  the  values  of  <pv  which  one  would  obtain  at  point  K if  one  solves 
the  original  difference  equation  (by  the  marching  procedure) . A 
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curve  drawn  through  the  points  <J>  will,  in  general,  have  a break  at 
the  point  K-l.  The  starting  approximation  cannot  be  expected  to 
have  the  correct  value  of  for  this  value  is  not  known  in 

advance.  The  contribution  of  the  homogeneous  particular  solutions 
to  the  starting  solution  will  therefore  be  finite  at  the  point  K-l, 
while  it  is  zero  at  point  K.  One  sees  that  homogeneous  particular 
solutions  of  short  wave  length  will  play  an  important  role  in  the 
hyperbolic  problems.  In  elliptic  equations,  a corresponding  jump 
tends  to  zero  as  the  mesh  is  refined. 
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SECTION  V 


STABILITY 


The  stability  of  the  procedure  for  a certain  choice  of  the 
parameters  ot^,  a 2»  , and  n/N  can  be  studied  by  assuming 

arbitrary  initial  conditions  for  <J>  which  satisfy  the  boundary 
conditions  of  the  homogeneous  problem  and  by  carrying  the  iteration 
out.  The  system  of  equations  for  the  is  actually  a three- 

point  recurrence  relation  with  constant  coefficients.  It  can 
be  efficiently  solved  by  the  method  of  cyclic  reduction. 

(Because  of  this  method  we  have  chosen  K = 32  = 25) . Such  compu- 
tations have  been  carried  out  for  a number  of  cases.  The  <t>j^ 
have  been  plotted  for  different  values  of  i.  A sample  of  these 
plots  is  shown  in  Figures  1 and  2.  The  initial  conditions  chosen 
are  ct>k  = 0 for  all  points  except  either  the  point  k = 1 or 
k = K-l.  Such  initial  conditions  could  be  represented  by  a 
linear  combination  of  the  eigensolutions  of  the  problem  (with 
none  of  the  coefficients  being  zero) . Therefore,  it  cannot  happen 
that  an  eigensolution  which  might  cause  instability  (that  is  for 
which  the  eigenvalue  has  an  absolute  value  larger  than  one)  is 
inadvertently  excluded  from  the  initial  condition.  In  spite  of 
this  fact,  one  finds  practically  no  similarity  in  the  convergence 
behavior  for  these  two  kinds  of  initial  conditions. 

A more  detailed  insight  is  obtained  by  computing  the  eigen- 
functions and  eigenvalues  which  determine  the  iteration  process. 

Let 


. fc)  i 

A * % * > 


km  'j 


K-t 


The  u.  's  are  combined  into  a vector  u.  Then  one  obtains  from  Eauation 
k 

(7)  the  following  matrix  eigenvalue  problem. 


[ A + $]  £ - £ [A  + C] U - D 


A,  B and  C are  matrices  of  dimension  K-l  by  K-l.  Our  examples 
have  been  computed  for  K = 32. 


with  a = 1+  0^ 

b = - (1-a  1)  -(1+  o2) 
c = (1  -ct±) 

C = - 4 sin2(ir^  ) I31 


where  is  the  31  by  31  identity  matrix 


In  our  computations  the  15*"^  point  has  been  taken  as  the  shock 
point.  The  eigenfunctions  and  eigenvalues  are  computed  using  the 
Q-Z  Algorithim  Ref.  (4,  5)*. 

The  case  | u^ | = .5;  n/N  = 0.05;  = 0;  a ^ = 1 has  been  taken 

as  the  starting  point  of  the  explorations. 

For  K = 32  one  has  31  eigenvalues  and  eigenfunctions,  most 
of  them  occur  in  conjugate  complex  pairs.  Figures  3 and  4 show  for 
the  elliptic  and  the  hyperbolic  problem  the  eigenvalues  which  lie 
in  the  first  quadrant  of  the  complex  p plane.  They  lie  rather 
close  to  each  other.  Accordingly,  one  will  not  have  an  outspoken 
dominance  of  a single  eigenvalue  (or  rather,  in  this  case,  of  a 
pair  of  conjugate  complex  eigenvalues > even  after  a fairly  large 
number  of  iterations.  Figures  6 through  17  show  some  of 
the  eigenfunctions  for  the  elliptic  and  for  the  hyperbolic  cases. 
They  are  ordered  according  to  the  absolute  value  of  p . The 
waviness  of  the  eigenfunctions  increases  as  the  absolute  value 
of  p decreases.  It  is  somewhat  disturbing  that  at  the  left  end  of 
interval  the  eigenfunctions  are  extremely  small.  To  represent  an 
initial  condition  in  which  all  values  of  < are  zero  except  for 
(j>^,  requires  a linear  combination  of  the  eigenfunctions  with 
very  large  coefficients.  In  contrast  if  <}>„  .=  1 and  all  others 
are  zero  the  coefficients  of  this  linear  combination  will  be  of 
moderate  size.  This  explains  the  difference  of  the  results  in 
Figures  1 and  2. 

In  Figures  3 and  4 some  values  of  p which  hold  for  an 
infinite  interval  are  included.  To  obtain  some  similarity  of 
the  problem  for  an  infinite  interval  with  the  eigenvalue  problem 
for  a finite  region,  we  assumed  that  one  has  the  same  interval 
but  with  periodicity  conditions  for  k = 0 and  k = 32.  This 
then  selects  32  discrete  eigenvalues  from  the  continuous  set 
given  by  Equation  (9) . The  values  of  y in  Equation  (9)  are 

*The  algorithm  is  part  of  the  EISPACK  library,  but  it  is  not 
yet  contained  in  the  second  edition  of  the  EISPACK  Guide. 
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found  from  the  periodicity  requirement 

hence 


where 


-KA  < 


* J k 

m-  v<  v-  K/i- 


The  eigenvalues  for  periodicity  conditions  lie  in  the  vicinity  of 
those  obtained  with  <J>  = 0 at  both  ends  of  the  region,  but  not  in  a 
close  vicinity.  The  discrepancies  are  particularly  large  for  the 
longest  and  for  the  shortest  wave  lengths.  To  recognize  the 
influence  of  the  length  of  the  region  under  consideration,  we  have 
carried  out  the  computations  with  boundary  conditions  <j)  = 0 also 
for  K = 16  (Fig.  5) . All  eigenvalues  so  obtained  are  close  to  some 
of  the  eigenvalues  of  the  problem  for  K = 32.  One  sees  that  the 
stability  analysis  for  an  infinite  interval  gives  results  which 
are  different  from  those  for  a finite  interval  (even  though  32  is 
a rather  large  number  ).  In  the  present  example,  the  stability 
criteria  obtained  for  an  infinite  interval  are  conservative. 


Figure  18  shows  the  eigenvalues  for  an  elliptic,  a hyper- 
bolic operator,  and  a problem  with  a shock.  The  number  of  eigen- 
values is,  of  course,  the  same  for  the  three  cases.  In  this  par- 
ticular case  one  finds  the  remarkable  result  that  the  eigenvalues 
for  the  shock  are  close  either  to  elliptic  or  hyperbolic  eigen- 
values. The  same  trend,  but  less  pronounced,  is  found  in  Figure  18 
for  n/N  =0.  In  any  case,  one  observes  that  the  presence  of  the 
shock  has  no  adverse  effect  on  the  stability  of  the  procedure. 

Figures  19,  20  and  21  give  some  indication  of  what  happens 
if  the  parameters  of  the  problem  are  changed.  This  part  of 
the  study  is  by  no  means  exhaustive.  One  sees  that  a decrease 
in  n/N  and  of  u^  has  a destabilizing  effect,  an  increase  in 
(*2  makes  the  problem  less  stable  in  the  elliptic  case  and 
more  stable  in  the  hyperbolic  case. 


— 
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These  graphs  can  be  used  to  discuss  the  influence  of 

over  or  under  relaxation  between  the  relaxation  steps.  Let  u, 

th  ^ 

be  the  k eigenvector  and  let  the  difference  between  the  values 

of  <j>  in  a certain  iteration  step  and  in  the  final  solution  be 


expressed  as 


u^.  Then  one  obtains  in  the  next  iteration 


step  for  this  difference  ^ C Pv  u,  . The  change  of  this 


— '-'-'jr  — v — w**_*.~  ^ ]<^  f ]<i.  k*  v/uauyc  wj.  Lino 

k=1  K-l 

difference  from  one  step  to  the  next  is  given  by  C^d-P^u^. 

k=l 

If  one  uses  a relaxation  factor  |3  then  this  change  is  applied  only 
partially,  that  is,  one  takes  as  a result  of  the  next  iteration 
step  i<_,  /(„,  ^ t<_, 

4*/  ' 

For  a relaxation  with  a relaxation  factor  B,  one  therefore  has 
as  effective  eigenvalues 

One  will  try  to  choose  3 in  such  a manner  that  m^x  | | is  as 
small  as  possible;  3 is  of  course  real.  Since  ,for  the  critical 
eigenvalues^  p^  - 1 is  nearly  imaginary,  one  can  attain  only  very 
little  improvement  by  the  choice  of  a relaxation  factor  3. 


SECTION  VI 


I 


THE  "SECOND"  METHOD  OF  MARTIN  AND  LOMAX 


Martin  and  Lomax  recommend  an  extrapolation  procedure  for 

the  results  of  subsequent  iteration  steps  based  on  the  Aitken- 

Shanks  formula.  This  formula  is  correct  if  the  d.'  ‘ferences  between 

the  results  of  the  iteration  steps  behave  like  a geometric  series. 

To  be  specific,  let  the  results  of  3 subsequent  iterations  be 
12  3 

given  by  <f>  , , <f>  , , and  <f>  and  assume  that  the  differences 

.K  K.  K. 

2 1 3 2 

(<j>  k “ <f>  k^  ' (<1>  k “ ^ k^  also  subsequent  differences)  form 

a geometric  series.  Then  one  obtains 


0 *yif 

oo 


(10) 


This  is  the  Ai tken-Shanks  formula.  Accordingly,  starting  with  an 
initial  value  , one  carries  out  two  iteration  steps  and  then 
uses  this  formula  to  extrapolate.  This  gives  a new  initial  con- 
dition which  is  then  treated  in  the  same  manner.  The  differences 
in  subsequent  iteration  steps  can  be  expressed  by  the  eigenfunctions 

u,  . The  individual 
k 

terms  of  the  sums  are  given  by  geometric  series.  If  there  were  a 
predominant  eigenvalue,  then  one  would  obtain  a behavior  of  the 
differences  close  to  that  of  a geometric  series  after  a number  of 
iteration  steps,  and  one  would  expect  this  formula  to  be  effective. 
For  the  distribution  of  eigenvalues  obtained  under  the  present 
circumstances,  this  is  unlikely  because  the  eigenvalues  are  so 
closely  spaced.  An  attempt  to  apply  the  procedure  in  the  present 
problem  was  without  success.  The  convergence  with  the  use  of  the 
Ai  tken-Shanks  formula  \yorse  than  that  of  direct  iterations. 


and  eigenvalues:  EC^  u^,  EC^  ^ u^ , EC^.  p* 
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It  is 


In  addition,  we  have  used  an  improved  version, 
based  on  the  observation  that,  in  the  present  problem,  the 
eigenvalues  occur  in  complex  conjugate  pairs.  Under  favorable 
circumstances,  two  eigenvalues  will  be  dominant.  The  general 
version  of  the  Aitken-Shanks  formula  holds  if  the  differences 

^k  - ^k ' ^k  “ ^k'  ^k  ” ^k  ***  are  sums  of  the  corresponding 

terms  of  a finite  number  of  geometric  series.  This  observation 
can  be  found  in  the  Section  "Heuristic  Motivation  of  the  Transform" 
of  Ref.  5 . 

Following  Shanks,  we  denote  A(n=0,l,2,...)  a sequence 
of  numbers  or  functions  and  set 


AA 


i 


Ai  + 1 


(11) 


The  starting  point  is  Equation  2 of  Reference  5.  For  the 
present  purpose,  we  set  n = k.  The  extrapolated  value  (that  is 
the  value  which  one  obtains  by  applying  the  Aitken-Shanks  formula) 
is  then  given  by 


Bk,k 


VD2 


where 


(12) 


(13) 


(14) 


and  M is  the  matrix  of  dimension  k by  k+1  whose  elements  are  given 
by 


M . . = AA . , . 

1/3  1+3 


i = 0,1,2... k-1 
j = 0 , 1 , 2 . . .k 


(15) 
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We  express  the  A^  by  the  first  term  of  the  sequence  A^  and  the 
differences  A A^.  One  has  for  the  row  vector  occurring  in 


E aq ' A1 ' A2 ' * * * Av  ^ ~~  E 1 / 1 r 1 ; • • • 1 ] t [ 0 , A1 -An , A^ -An , • • • A^  ~A n ] 


1 "O'  2 O' 


k 0J 


Hence,  from  Equation  (12), 


Bk,k  A0  + D3/D2 


where 


D3  * 


A2_A0 ' * * * Ak  A0 


M 


-) 


Let  us  evaluate  under  the  assumption  that 


(16) 


A A. 


= Vb  rl  , 

/ j n n 


£ = 0,  1,  2 . . . 


n=i. 


(17) 


The  r 's  are  the  ratio  of  consecutive  terms  in  the  individual 
n 

geometric  series.  Then  one  obtains 

*n-i  £ ■**-/  4 ^ %£ 


A -A 

m,  t> 


£ - Z 4,  ? 


i--o 


c *v=i  e=o 


*V=! 


? * 


(18) 


Therefore , 


[0,  A^  Aq,  A2  Aq , . . . , A^-Aq ] 
k b 

2 1?^  { [1,1,1,  ...  ,11  - [1,  rn,  rn  ,...,rnk]| 

n=l 

This  is  substituted  into  Equation  (16) . Using  Equation  (14) , one 
obtains 

k b 


n 


Bkk=  A0  + 


1-r . 


n 


V°2 


n=l 


(19) 
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where 


£ — 
/~r*. 

-r  ^ r*v 

T 4+ 

~ /- 

• Z ^ ^ 

' - r*v 

£ 

/4C  •• 

• r+„  --J 

2 * 

•h 

£ *5, 

• 

£ 4 * ’ 

• 

• 

# 

<t-> 

£ Cc 

• 

7 4 
*->  ***'*' 

• 

Z*/  ■ 

• ■ ^t4. 

The  rows  of  the  k+1  by  k+1  determinant  can  be  expressed  as  linear 
combinations  of  only  k vectors;  which  appear  as  the  rows  of  the 
third  matrix  in  the  following  equation: 


Therefore,  D„  = 0 

4 


and  k 

Bkk  * A0  + Yj  v(i-v 

n=l 
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With  Equations  (17)  and  (11)  one  then  obtains 


To  take  the  fact  into  account  that  the  dominant  eigenvalues 
of  our  problem  occur  in  conjugate  complex  pairs  we  have  used  the 
Aitken  Shanks  formula  for  k = 2.  Then  one  deals  with  3 by  3 
determinants,  and  the  difference  between  consecutive  terms  of  the 
sequences  are  approximately  the  sum  of  the  respective  terms  of 
two  geometric  series.  One  has  specifically 


4 


/ 

/ 

/ 

O 

aA0 

(aAc  vaA,) 

& A0 

*4, 

A At_ 

j A, 

a 

4 

(22) 


a,  ^ 

(a  A,  -a  -/aAe)aAs  - a A,a  AJ  + /aAoi j 4^-aA,) 


(23) 

* - (*\*\-*\*  A±-)  + aA,)(aA^  aA2  - jA*) 


(24) 


Under  certain  conditions,  the  determinant  D2 , Equation  (14)  may 
vanish.  Assuming  that  A A^  has  the  form  (16) , one  has 


/ 

/ 

/ • • . 

/ 

Z A r 

Z-6-  r1  • • ■ 

«• 

■ ■ ■ 

Z<L 

Z < 

■ ■ 

■ Xiv*" 

Z A-  r*’’'  Z A-  r* 

J1  . 

_,  . z4- 

• 
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The  matrix  encountered  here  can  be  written  as  the  product  of  three 
k+1  by  k+1  matrices 


It  is  evident  that  D-,  will  vanish  if  two  or  more  values  r coincide, 

c n 

or  if  some  of  the  r ' s are  zero.  Incidentally,  the  first  and  third 
of  these  matrices  are  closely  related  to  Vandermond ' s determinant. 

On  this  basis,  one  obtains: 

* '£*••  h)lh- Vlt(r3- Vfy 

~ £ -i  )(rk,  ' £-3  ) ’ ' ’ f\~  7i)U  x 

* Ufr'  ][Cr3~\  )(r3  - -$//£  - V(rf  - rjfy 

(25) 

If  one  or  more  of  the  r^'s  are  zero,  then  the  sum  of  the  corresponding 
geometric  series  gives  infinity.  The  extrapolation  given  by 
Equation  (12)  or  (24)  is  therefore  correct,  although  one  will 
question  whether  the  formula  should  be  applied  under  these 
circumstances . 

The  following  discussions  are  restricted  to  the  case  k = 2. 

For  = r^  / 1 one  has  just  one  geometric  series  rather  than  two. 

In  this  case,  one  obtains  0/0  from  Equation  (24)  . One  must  go 
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back  to  the  formula  for  k = 1.  This  case  is  characterized  by 
the  conditions 

A A0  A A2  - A Ax2  = 0 


A A±  A Aq  - A A2  = 0 


The  result  is  then  given  by 

A0  + A A0/<1  - A Aj/A  A0) 


AA.  AA,  - AA  / 0 


AA^  AA^  - AA2  = 0 


then  one  of  the  values  of  r is  zero  and  one  has 


AA0  = bi  + b2 


AAX  = 


AA  = 
2 


b2r2 

b2r2 


aa3  = 


b2r2 


In  this  case.  Equation  (24)  is  applicable,  but  one  must 
decide  whether  vanishing  of  one  of  the  r's  is  compatible  with  the 
nature  of  the  problem. 


A AQ  AA2  - AAX  = 0 


AAX  AA3  - AA2  / 0 


then  one  of  the  r's  is  infinite.  One  has 
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AA0  = b2 


AA1  = b2r2 


AA2  = t»2r2 
AA3  = b2r2  + b1 


(aa4  = ») 


This  case  must,  of  course,  be  excluded  although  Equation  (24) 
will  assign  a limiting  value  to  the  sequence. 

The  case 


AArt  = c. 


AA,  = 


C1  rl  + C2 


AA~  = 


C1  rl  + 2 c2  rl 


AA.,  = c.  r,  + 3 c~  r 


'1  "1 


2 L1 


arises  by  a limiting  process  where  r2  -*■  r3  and  b2  -*•  °°,  b.^  00 . 

Here  Equation  (24)  is  applicable  without  modifications. 

We  have  experimented  with  Equation  (24)  as  a means  of 
accelerating  the  convergence  of  the  iterative  procedure  whose 
results  are  shown  in  Figures  1 and  2.  In  this  case,  one  generates  five 
consecutive  values  of  <j>  ^ and  then  uses  (24)  to  extrapolate. 

The  results  give  initial  values  for  a next  step  of  the  same  kind. 

It  was  found  that  the  results  are  somewhat  better  than  those  obtained 
with  Equation  (10) , but  there  is  no  noticeable  improvement  of  the 
convergence  in  comparison  with  straight  iterations. 

One  will  remember  that  the  extrapolation  procedure  is  a 
non-linear  process.  It  is  therefore  difficult  to  generalize 
results  obtained  by  examples. 
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The  practical  aspects  of  the  Martin-Lomax  iterations  have  been 
quite  thoroughly  explored  in  the  work  of  these  authors.  To  this 
facet  of  the  work,  the  present  discussions  cannot  claim  to  make 
a contribution.  We  hope  that  they  are  useful  as  background 
information  in  cases  where  one  encounters  phenomena  that  are 
difficult  to  explain.  In  this  regard,  the  data  about  eigenfunctions 
and  eigenvalues  can  be  rather  revealing.  The  observation  that  the 
presence  of  a shockpoint  does  not  lead  to  a significant  change  of 
the  convergence  properties  is,  of  course,  reassuring.  The 
discussion  (carried  out  in  Section  I)  of  the  relation  between  the 
boundary  conditions  for  a hyperbolic  problem  and  those  used  in  the 
present  iteration  scheme  is  important  for  the  understanding  of  the 
working  of  the  procedure.  According  to  the  present  examples,  it 
seems  to  be  rather  difficult  to  find  values  of  the  free  parameters 
and  a ^ which  guarantee  convergence  under  all  circumstances. 
Martin's  numerical  experiments  show  that  these  difficulties  are  less 
pronounced  if  the  procedure  is  applied  to  computations  of  an  actual 
flow  field.  The  fact  that  the  eigenvalues  p are  clustered  closely 
around  1 (at  least  in  the  present  examples)  is  disappointing,  it  may 
mean  that  convergence  is  somewhat  precarious;  this  may  create 
difficulties  if  one  wants  to  make  procedures  of  this  kind  fully 
automatic. 
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ial  conditions: 


ial  conditions 


Figure  5. 


Eigenvalues  in  the  complex  p plane 
|u^|  = .5,  n/N  = 0.05,  = 0,  = 

full  small  circles  K = 32 
open  circles  K = 16 


1 


Figure  7.  Real  and  imaginary  part  of  eigenvector  for 
n/N  = .05,  = 0,  = 1,  hyperbolic  case 

Pr  = 0.88120  + 0.22250  i. 


Figure  8.  Real  and  imaginary  part  of  eigenvector 
n/N  = .05,  = 0,  a ^ = 1,  hyperbolic  < 

PQ  = 0.79986  + 0.17233  i. 


Figure  9.  Real  and  imaginary  part  of  eigenvecto 
n/N  = .05,  = 0,  = 1/  hyperbolic 

p, , = 0.72738  + 0.06853  i. 


Figure  11.  Eigenvector  for  | | = .5,  n/N  - .05,  a.  - 0, 

a = 1,  hyperbolic  case.  p=  0.64256  (real),  real. 


hypervolic  case,  p = 0.55762  (real)  , u,  real 


Hip 


<N 

m 


CM 

co 


.68359  + 0.16305 


Figure  15.  Real  and  imaginary  part  of  eigenvector  for 
n/N  = .05,  = 0,  = 1,  elliptic  case 

pQ  = 0.65600  + 0.08698  i. 


p = 1 

Figure  18.  F.igenvalues  for  finite  interval, 

hyperbolic  problem  empty  circles 

elliptic  problem  empty  squares 

problem  with  shock  small  full  circles 

| u | = .5,  cu  = 0 , a_  = 1,  upper  graph  n/N 

lower  graph  n/N 


p = 1 


hyperbolic  case 


complex  p-plane 


□ 

□ 

a A 

□ o 


P = 1 


p = l 


elliptic  case 


complex  p-plane 


complex  p-plane 


complex  p-plane 


* - ■-  ■ - ■ ■■  a 

Figure  21.  Eigenvalues  in  the 

P = 1 

complex  p plane 

t 

n/N  = 0.05,  | u£  | = 

.5,  = 0, 

'f  f 

empty  circles 

a2  = 2 

full  small  circles 

a2  = 1.5 

empty  squares 

a2  = 1 
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